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Abstract 

We solve the nonequilibrium dynamics of a 3 + 1 dimensional 
theory with Dirac fermions coupled to scalars via a chirally invariant 
Yukawa interaction. The results are obtained from a systematic 
coupling expansion of the 2PI effective action to lowest non-trivial 
order, which includes scattering as well as memory and off-shell effects. 
The dynamics is solved numerically without further approximation, 
for different far-from-equilibrium initial conditions. The late-time 
behavior is demonstrated to be insensitive to the details of the 
initial conditions and to be uniquely determined by the initial energy 
density. Moreover, we show that at late time the system is very well 
characterized by a thermal ensemble. In particular, we are able to 
observe the emergence of Fermi-Dirac and Bose-Einstein distributions 
from the nonequilibrium dynamics. 
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1 Introduction 



The abundance of experimental data on matter in extreme conditions 
from relativistic heavy-ion colhsion experiments, as well as applications 
in astrophysics and cosmology has lead to a strong increase of interest 
in the dynamics of quantum fields out of equilibrium. Experimental 
indications of thermalization in collision experiments and the justification of 
current predictions based on equilibrium thermodynamics, local equilibrium 
or (non-)linear response pose major open questions for our theoretical 
understanding One way to resolve these questions is to try to 

understand quantitatively the far-from-equilibrium dynamics of quantum 
fields, without relying on the assumption of small departures from 
equilibrium, or on a possible separation of scales that forms the basis 
of effective kinetic descriptions 0. In contrast to close-to-equilibrium 
approaches, the quantum-statistical fluctuations of the fields are not assumed 
to be described by a thermal ensemble. Moreover, one goes beyond the range 
of applicability of the usual gradient expansion and dilute-gas approximation. 

In contrast to thermal equilibrium, which keeps no information about 
the past, nonequilibrium dynamics poses an initial value problem: time- 
translation invariance is explicitly broken by the presence of the initial 
time, where the system has been prepared. The question of thermalization 
investigates how the system effectively looses the dependence on the details of 
the initial condition, and becomes approximately time-translation invariant 
at late times. According to basic principles of equilibrium thermodynamics, 
the thermal solution is universal in the sense that it is independent of the 
details of the initial condition and is uniquely determined by the values of 
the (conserved) energy density and of possible conserved charges.^ There are 
two distinct classes of universal behavior, corresponding to Bose-Einstein and 
Fermi-Dirac statistics respectively. 

The description of the effective loss of initial conditions and subsequent 
approach to thermal equilibrium in quantum field theory requires calculations 
beyond so-called "Gaussian" (leading-order large-A^, Hartree or mean-field 
type) approximations j^Elini- Similar to the free-field theory limit, these 
approximations typically exhibit an infinite number of additional conserved 
quantities which are not present in the underlying interacting theory [71 |H] . 

^Here we consider closed systems without coupling to a heat bath or external fields, 
which could provide sources or sinks of energy. 
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These spurious constants of motion constrain the time evolution and lead to 
a non-universal late-time behavior P]. It has recently been demonstrated in 
the context of scalar field theories, that the approach to quantum thermal 
equilibrium can be described by going beyond these approximations IH] . 
In particular, this has been achieved by using a systematic coupling expansion 
[TUj or a expansion to next-to- leading-order jHl E] of the two-particle 
irreducible generating functional for Green's functions, the so-called 2PI 
effective action O ESI Ell CSl • In this context, the corresponding equations 
of motion have been shown to lead to a universal late time behavior, in 
the sense mentioned above, without to have recourse to any kind of coarse- 
graining. 

In this work, we study the far-from-equilibrium time evolution of 
relativistic fermionic fields and their subsequent approach to thermal 
equilibrium. Nonequilibrium behavior of fermionic fields has been previously 
addressed in several scenarios IHl UH CZl UH Ull; including the full 
dynamical problem with fermions coupled to inhomogeneous classical 
bosonic fields [201 ■ Here we go beyond these approximations and compute the 
nonequilibrium evolution in a 3-1-1 dimensional quantum field theory of Dirac 
fermions coupled to scalar s in a chirally invariant way. As a consequence, 
we are able to study the approach to quantum thermal equilibrium. For 
the considered nonequilibrium initial conditions we find that the late-time 
behavior is universal and characterized by Fermi-Dirac and Bose-Einstein 
distributions, respectively. The results are obtained from a systematic 
coupling expansion of the 2PI effective action to lowest nontrivial (two- 
loop) order, which includes scattering as well as memory and off-shell 
effects. The nonequilibrium dynamics is solved numerically without further 
approximations. We emphasize that, given the limitations of a weak-coupling 
expansion, this is a first-principle calculation with no other input than the 
dynamics dictated by the considered quantum field theory for given initial 
conditions. 

In Sect. we review the 2PI effective action for fermions, which we use 
to derive exact time evolution equations for the spectral function and the 
statistical two-point function in Sect. El We discuss the Lorentz structure of 
our equations in Sect. Eland exploit some symmetries in Sect. El In Sect. 
we specify to a chiral quark model for which we solve the nonequilibrium 
dynamics. The initial conditions are discussed in Sect, d and some details 
concerning the numerical implementation are given in Sect.jHl The numerical 
results are presented and discussed in Sect. (HI and Sect. |TU1 We attach two 
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appendices discussing in detail the quasiparticle picture we used to interprete 
some of our results. 



2 2PI effective action for fermions 

We consider first a purely fermionic quantum field theory with classical action 

S = J dV(v^i(x)[«^ -m/]^i(x) + \/(V^,^)) (2.1) 

for i = 1, . . . , Nf "flavors" of Dirac fermions ipi, a mass parameter rrif and an 
interaction term Vltpjtp) to be specified below. Here ^ = 'J^df^, with Dirac 
matrices 7^ (/x = 0, . . . , 3). Summation over repeated indices and contraction 
in Dirac space is implied. All correlation functions of the quantum theory can 
be obtained from the corresponding two-particle-irreducible (2PI) effective 
action F. For the relevant case of a vanishing fermionic "background" field 
the 2PI effective action can be written as [13 

T[D] = -iTilnD-^ - iTi Dq^D + T2[D] + const. (2.2) 

The exact expression for the functional F2 [D] contains all 2PI diagrams with 
vertices described by V{ip,ip) and propagator lines associated to the full 
connected two-point function D. In coordinate space the trace Tr includes 
an integration over a closed time path C along the real axis :21J, as well as 
integration over spatial coordinates and summation over flavor and Dirac 
indices. The free inverse propagator is given by 

iDo}j{x, y) = {i^ - ruf) 5^(x - y) 5ij . (2.3) 

The equation of motion for D in absence of external sources is obtained by 
extremizing the effective action jTS] 

6T\D] , , 

^ ^ -0. (2.4) 



bDi^{x,y) 

According to ()2.2|) one can write ()2.4|) as an equation for the exact inverse 
propagator 

D-^{x, y) = D^}^{x, y) - S,,(x, y- D) , (2.5) 
with the proper self-energy 

T.,^{x,y-D)^-i (2.6) 
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Equation (j2.5|) can be rewritten in a form, which is suitable for initial value 
problems by convoluting it with D. One obtains the following time evolution 
equation for the propagator: 

[i^^ - mf)Dij{x, y) -i z; D)Dkj{z, y) = i6c{x - y)6ij , (2.7) 

J z 

where we employed the shorthand notation = dz^ J d^z. Note that to 
keep the notation clear, we omit Dirac indices and reserve the Latin indices 
k, ... to denote flavor. 

3 Exact evolution equations for the spectral and 
statistical components of the two-point function 

To simplify physical interpretation we rewrite ()2.7p in terms of equivalent 
equations for the spectral function, which contains the information about the 
spectrum of the theory, and the statistical two-point function. The latter will, 
in particular, provide an effective description of occupation numbers. For 
this we write for the time-ordered two-point function Dij{x,y) and proper 
self-energy T,ij{x, y; D):^ 

D,jix,y) = ecix'-y')D>^ix,y)-eciy'-x')D<^{x,y), (3.1) 
j:i,{x,y;D) = ec(xO-2/°)S>(a;,y)-ec(/-a;0)S<(a;,y). (3.2) 

Note that for convenience we omit the explicit D-dependence in the notation 
for Sj^'^. Inserting the above decompositions into the evolution equation 
(12. 7p . we obtain evolution equations for the functions Dfj{x, y) and D^{x, y) 
as well as the identity 

7° {Dlj{x,y) + D<{x,y)) \^o=yO = 6{x-y)6,j, (3.3) 
^If there is a local contribution to the proper self-energy, we write 

E,,(x, y; D) - -i D) S^ix - y) S,, + y; D) , 

and the decomposition H3.2() is taken for y; D). In this case the local 

contribution gives rise to an effective space-time dependent fermion mass term mf + 

J.(l0Cal)(2;;£)). 
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which corresponds to the anticommutation relation for fermionic field 
operators. For later use we also note the hermiticity property 

{D>iy,x)y = ^'D>^ix,yh\ (3.4) 

and equivalently for Df-{x,y). Note that here the hermitean conjugation 
denotes complex conjugation and taking the transpose in Dirac space only. 

Similar to the discussion for scalar fields in Refs. 122] we introduce the 
spectral function pij{x,y) and the statistical propagator Fij{x,y) defined as^ 

p,,{x,y) = t(^D>^{x,y) + D<^{x,y)y (3.5) 

F,^{x,y) = ^(^D>^{x,y)-D<^{x,y)y (3.6) 
The corresponding components of the self-energy are given by^ 

A.,(x,i/) = ^(S>.(x,y) + S<.(a;,y)), (3.7) 

C^,{x,y) = i(s>(x,y)-S<.(x,y)). (3.8) 

With ()3.4j) the two-point functions have the properties 

iPj^iy,x))^ = -7%(x,y)7°, (3.9) 
{Fj,{y,x))^ = -f^F,,{x,yh\ (3.10) 

and equivalently for Aij{x,y) and Cij{x,y). 

Using the above notations the evolution equation ()2.7|) written for Pij{x, y) 
and Fij{x, y) are given by 



('^^x-^f)piji^^y) = / dzA,k{x,z)pkj{z,y), (3.11) 

JyO 

{i^^-mf)Fij{x,y) = / dz Aik{x, z)Fkj{z,y) 

Jo 

- / dzCik{x,z)pkj{z,y) , (3.12) 
Jo 



^Equivalently, one can decompose 



^Besides the dynamical field degrees of freedom Dij we introduce quantities which are 
functions of these fields. These fimctions are denoted by either boldface or Greek letters. 
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where we have taken the initial time to be zero, and J^o dz = J^o / d^z. 
For known self-energies the equations (j3.1ip and (|3.12|) are exact. We note 
that the form of their RHS is identical to the one for scalar fields [HI 1^ . 
To solve the evolution equations one has to specify initial conditions for 
the two-point functions, which is equivalent to specifying a Gaussian initial 
density matrix^. We note that the fermion anticommutation relation or (|3.3|) 
uniquely specifies the initial condition for the spectral function: 

'y"pij{x,y)\^o=yO =i5{x-y)5ij. (3.13) 

Suitable nonequilibrium initial conditions for the statistical propagator 
Fij{x,y) will be discussed below. 

4 Lorentz decomposition 

It is very useful to decompose the fields Pij{x, y) and Fij{x, y) into terms that 
have definite transformation properties under Lorentz transformation. We 
will see below that, depending on the symmetry properties of the initial state 
and interaction, a number of these terms remain identically zero under the 
exact time evolution, which can dramatically simplify the analysis. Using a 
standard basis and suppressing flavor indices we write 

P = Ps + h5pp + If^Pv + If^lbPA + ^(^f^upT ' (4-1) 

where cr^jy = |[7^,7i/] and 75 = i'j'^'y^'y^'y^ . For given flavor indices the 16 
(pseudo-) scalar, (pseudo-) vector and tensor components 

Ps = trp, 

pp = -itij^p, 

p'v = tr7V, (4.2) 

Pa = tr757''P, 
p^^ = tra^V, 

""We emphasize that a Gaussian initial density matrix only restricts the initial conditions 
or the "experimental setup" and represents no approximation for the time evolution. More 
general initial conditions can be discussed using additional source terms in defining the 
generating functional for Green's functions. 
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are complex two-point functions. Here we have defined tr = ^tr where the 
trace acts in Dirac space. Equivalently, there are 16 complex components for 
Fij, Aij and for given flavor indices Using ()3.9p and ()3.10p . one sees 
that they obey 

p!?(^,y) = -(p$?(y,^))* , F^Pix,y) = {F^Piy.x))\ (4.3) 

where F = {S, P,V, A,T}. Inserting the above decomposition into the 
evolution equations ()3.1H) and ()3.12|) one obtains the respective equations 
for the various components displayed in Eq. ()4.2j) . 

For a more detailed discussion, we first consider the LHS of the 
evolution equations 1)3.111) and ()3.12|) . In fact, the approximation of a 
vanishing RHS corresponds to the standard mean-field or Hartree-type 
approaches frequently discussed in the literature jUEj- However, to discuss 
thermalization we have to go beyond such a "Gaussian" approximation: it 
is crucial to include direct scattering which is described by the nonvanishing 
contributions from the RHS of the evolution equations. Starting with the 
LHS of ()3.11|) one finds, omitting flavor indices (see also Ref. 



tr [{i^ 


- mf)p\ 


= (id^Pv) -^fPs , 


-i tr [75 (i^ 


- mf)p] 


= -i{id^p^^) -nifpp , 


tr [7''(i^ 


- ^f)p] 


= (id'^ps) + i (id^pT^) -nifp'^ , 


tr[757'^(«^ 


- mf)p] 


= t (td^pp) + i e^""'' {td,pT,-y5) -mfp>X, 




- ^f)p] 


= -I {id^p\, - id'^p^y) + e^"^' {id,pA,5) -rrifp^^ 



(4.4) 



The corresponding expressions for the LHS of the evolution equation (|3.12p 
for F follow from (|4.4j) with the replacement p ^ F. Considering now the 
various component ()4.2|) of the integrand on the RHS of Eq. ()3.1H) . we find 

tr [A p] = Asps- Ap pp + py,^ - A^ p^,/. 

+iA^>T,,. , (4.5) 
-i tr [75 A p] = As Pp + Ap ps - i A^ pA,^, + i A^ py,^ 

+ ^e^'^^'AT,,uPT,,5 , (4.6) 
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+zAf py,y + ie'^"^' (Aa,. PTaS + At,.^ Pa,^) , (4.7) 
tr [757^^ A p] = Asp'X + K Ps - i^p Py + Pp + ^^a,u p7 

+z A^^ PA,u + ^e^-^^' ( Av^,, pT,^s + At,.^ py,^) , (4.8) 

tr y^'A p] = As pf^r + A^'^ PS - ^ ( Ap Pt,^5 + Ar,-,^ Pp) 

-2 (A^ p^^;, - A^ p^) + e^'^^'^ (Ay,^ pA,s - Aa,, pv,s) 

+z (a:^ p:^ - a:^ p^;) + z ( a^^ p^,/ - a^^ p^,/) . (4.9) 

With the above expressions one obtains the evolution equations for the 
various Lorentz components in a straightforward way using ()3.1ip . We note 
that the convolutions appearing on the RHS of the evolution equation (j3.12|) 
for F are of the same form than those computed above for p. The respective 
RHS can be read off Eqs. ()4.5|) - ()4.9p by replacing p ^ F for the first term 
and A — >■ C for the second term under the integrals of Eq. ()3.12|) . We 
have now all the relevant building blocks to discuss the most general case 
of nonequilibrium fermionic fields. However, this is often not necessary in 
practice due to the presence of symmetries, which require certain components 
to vanish identically. 

5 Symmetries 

In the following, we will exploit symmetries of the action ()2.H) and of 
the initial conditions in order to simplify the fermionic evolution equations 
derived in the previous section. 

Spatial translation invariance and isotropy: We will consider spatially 
homogeneous and isotropic initial conditions. In this case it is convenient to 
work in Fourier space and we write 

p{x,y)^p{x',y';x-y) = J e'P<'-y^ p{x' ,y';p) , (5.1) 

and similarly for the other two-point functions. Moreover, isotropy implies a 
reduction of the number of independent two-point functions: e.g. the vector 
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components of the spectral function can be written as 

where p = |p| and v = p/p. 

Parity: The vector components py{x^ ^y^]p) and pv{x^ ,y^',p) are 
unchanged under a parity transformation, whereas the corresponding axial- 
vector components get a minus sign. Therefore, parity together with 
rotational invariance imply that 

p°(a;°,l/°;p) = PA(x°,i/°;p) = 0. (5.2) 

The same is true for the axial- vector components of F, A and C. Parity also 
imphes the pseudo-scalar components of the various two-point functions to 
vanish. 

C P -invariance: For instance, under combined charge conjugation and 
parity transformation the vector component of p transforms as 

p'y{x\y'-p) pl{y\x^:p). 
Pv{x^.y^;p) — ^ -Py(l/°,a;°;p) , 

and similarly for Ay and Ay. The F-components transform as 
F^{x\y'-p) ^ 

and similarly for Cy and Cy. Combining this with the hermiticity 
relations ()4.3|) . one obtains for these components that 



(5.3) 



Rep^(x°,y°;p) = Impy(x°,?/°;p) = 0, 

ReF^(x°,?/°;p) = ImFy(a;°,?/°;p) = , 

ReA°.(x°,?/°;p) = Im Ay(a;°, /;p) = , 

ReC^(x°,?/°;p) = Im Cy(a;°, |/°; p) = , 

for all times x° and y^ and all individual momentum modes. 

We stress that a nonequilibrium ensemble respecting a particular 
symmetry does not imply that the individual ensemble members exhibit the 
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same symmetry. For instance, a spatially homogeneous ensemble can be build 
out of inhomogeneous ensemble members and clearly includes the associated 
physics. The most convenient choice of an ensemble is mainly dictated by 
the physical problem to be investigated. Below we will study a "chiral quark 
model" with Dirac fermions coupled to scalars in a chirally invariant way. 
Since in this paper we will restrict the discussion of this model to the phase 
without spontaneous breaking of chiral symmetry, it is useful to exploit this 
symmetry as well. 

Chiral symmetry: The only components of the decomposition (j4.H) 
allowed by chiral symmetry are those which anticommute with 75. We 
therefore have 

ps{x',y';p) = pp{x\y'-p) = pi;r{x^y'-p) = , (5.4) 

and similarly for the corresponding components of F, A and C. In particular, 
chiral symmetry forbids a mass term for fermions and we have nif = 0. 



5.1 Equations of motion 

In conclusion, for the above symmetry properties we are left with only four 
independent propagators: the two spectral functions py and pv and the two 
corresponding statistical functions Fy and Fy. They are either purely real 
or imaginary and have definite symmetry properties under the exchange of 
their time arguments ^ y^ . These properties as well as the corresponding 
ones for the various components of the self-energy are summarized below: 



^0 AO 

Pv^ 

Pv, Ay 
Fy, Cy 



imaginary, symmetric; 
real, antisymmetric; 
imaginary, antisymmetric; 
real, symmetric. 



The exact evolution equations for the spectral functions read (cf. Eq. ()3.11|) ): 
. d 



dx^ 



Pvi^ ;P) =PPv{x ,y ;p) 



(5.5) 



y" 



A^x", z"; p) pI{z\ y'; p) - Ay (x^ z"; p) pv{z\ y^- p) 



^Wc note that the foUowing equations do not rely on the restrictions (|5.3I) imposed by 
CP-invariance: they have the very same form for the case that aU two-point functions are 
complex. 
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d 



(5.6) 



Alix', z'- p) pv{z', y'- p) - Ay{x', z'- p) p^z", y'; p) 



.0 „,0 



^0. ^\ J) ( ^0 „,0. 



Similarly, for the statistical two-point functions we obtain (cf. Eq. ()3.12|) ): 

d 



t^F^ix'',y'';p)=pFvix'',y^;p) 



(5.7) 



+ / dz 
'o 



A^ix', z'- p) F^{z', y'- p) - Ay{x', z'- p) Fy{z', y'; p) 



^0., 



vO „.o. 



y 



dz 



C'yix', z'- p) pI{z\ y'; p) - Cy (x", z"; p) py{z' , y'; p) 



^0. 



.0 „,0. 



d 



i^Fy{x\y''-p)=pF^{x\y''-p) 



(5. 



+ / dz 

'0 



Al{x\ z'- p) Fyiz', p) - Ayix', z'; p) F^z', y'; p) 



.0 „,0 



^0. ^\ /^O „,0. 



dz 



C'y{x\ z'- p) py{z\ y"; p) - Qy{x\ z'- p) p%{z\ y"; p) 



The above equations are employed below to calculate the nonequilibrium 
fermion dynamics in a chiral quark-meson model. 



6 Chiral quark— meson model 

As an application we consider a quantum field theory involving two fermion 
flavors ( "quarks" ) coupled in a chirally invariant way to a scalar a-field and 
a triplet of pseudoscalar "pions" vr" (a = 1,2,3). The classical action reads 

+ gtp[a + Z75r"vr"] - V {a^ + ti^)] , (6.1) 



where vr^ = tt^tt'* and where t"' denote the standard Pauli matrices. The 
above action is invariant under chiral SUl{2) x SUr{2) transformations. For 
a quartic scalar self-interaction ~ (a^ + tt°'tt°') , this model corresponds to 
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the well known linear cr-model [201; which has been extensively studied in 
thermal equilibrium in the literature using various approximations [HI]. For 
simplicity we consider here a purely quadratic scalar potential: 



which is sufficient to study thermalization in this model. We note that this 
theory has the same universal properties than the corresponding linear o- 
model. Extending our study to take into account quartic self-interaction is 
straightforward. It gives additional contributions to the scalar self-energies, 
Eqs. fj6.13|) and ()6.14j) below. These contributions can be found in Refs. PlITT] 
and we will point out the respective changes below. 

6.1 Equations of motion for the scalar field 

The 2PI effective action for nonequilibrium scalar fields has been extensively 
studied in the literature O UHl UDl ISl HI]- Here, we briefly recall the main 
features of the scalar sector and stress those aspects which are relevant 
for the present paper (for details see e.g. Refs. El)- For the model 
considered here, Eq. ()6.H1 . the 2PI effective action is a functional of fermionic 
propagators as well as scalar propagators.^ The scalar fields form an 0(4)- 
vector 0A(a^) = (o"(a^), 7r(x) ) and we denote the full scalar propagator by 
Qab^^iV) with A,B = 0, . . . , 3. The 2PI effective action ()2.2|) . augmented 
by the scalar sector, reads 



Similar to the fermionic case discussed above, the equation of motion for the 
scalar two-point function is obtained by minimizing the 2PI effective action 
with respect to G, and the proper self-energy is given by [TSl El ITTj 



^As emphasized in the previous section, we do not consider the possibihty of a broken 
symmetry in this paper. Therefore we restrict the discussion to a vanishing scalar field 
expectation value. 




(6.2) 




^Go,ABi^^ y) = -i^x + ml) 5cix - y) 5ab ■ (6.4) 




(6.5) 
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Under chiral transformations, the matrix Q -^TZQTZ\ where TZ is an 0(4) 
rotation. Without loss of generahty, because of chiral symmetry the effective 
action T[Q, D] and its functional derivatives can be evaluated for Q taken to 
be the unit matrix in 0(4)-space: 

QAB{x,y) = G4,{x,y)5AB ■ (6.6) 

The same holds for the corresponding self-energy ()6.5|) . Similarly, in the 
fermionic sector, because of chiral symmetry the most general fermion two- 
point function can be taken to be proportional to unity in flavor space: 

Dij{x,y) = D{x,y)5ij . (6.7) 

Similar to the discussion in Sect. El the scalar spectral and statistical two- 
point functions are defined as p| 123 

G^(x, y) = F^{x, y) - '-p^{x, y) [0c(x° - /) - 0c(/ - x°)] , (6.8) 

and equivalently for the spectral and statistical self-energies and S^. 
These are all real functions and F-like components are symmetric under the 
exchange of x and y, whereas the p-like components are antisymmetric. The 
equal-time commutation relation of two scalar field operators implies jHl 1221 

P<l>{x, y) Uo=yO = 0, d^op^{x, y) \^o=yO =5{x-y), (6.9) 

which uniquely specifies the initial conditions for the spectral function. Initial 
conditions for the statistical two-point function will be discussed below. 
Finally, the equations of motion for the scalar propagators read |^ \'2'2\: 

[dlo+p' + ml]p^ix',y';p) = - f dz' ^{x\ z'^-p) p^{z\y'-p) . 

JyO 

(6.10) 

[dl,+f + ml]F^ix',y^;p) = - f dz^ ^x' , z^ ; p) F^iz' , y^ ; p) 

Jo 

+ r dz'T.l{x\z'-p)p^{z\y'-p), 
Jo 

(6.11) 
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where we have assumed spatially homogeneous initial conditions and Fourier 
transformed with respect to spatial coordinates. For known self-energies 
these are the exact equations for the theory described by the classical action 
()6.H) together with ()6.2p . In the presence of scalar self- interactions, the 
self-energy receives in particular a local contribution —i'E^^"^^^\x) 6^{x — y), 
which amounts to a shift of the bare mass squared appearing in the above 
equations. In that case, the exact equations of motion have the same form 
as above, with the replacement ml ^ M\x) =ml + Ej°'"')(a;) 0111]. 

6.2 2PI coupling expansion 

We consider a systematic coupling expansion of the 2PI effective action ()6.3p 
to lowest non-trivial order, which includes scattering as well as memory and 
off-shell effects, such that thermalization can be studied. The first non- 
trivial order in a coupling expansion corresponds to the two-loop contribution 
depicted in Fig.[T](cf. also the discussion in Ref. [221 Using and (j6.7p . 
one can express the two-loop contribution to r2 directly in terms of and D. 
We obtain for the chirally symmetric theory: 

^ -^/^ ^ d'x d'y ti[D{x, y) D{y, x)] G^{x, y) , (6.12) 

where Nf = 2 is the number of fermion flavors and A^^ = 4 is the number of 
scalar components. From there, it is straightforward to compute the spectral 
and statistical components of the two-loop self-energies^. We obtain for the 
scalar self-energies 

mx\y'-,p) = -^g'Nf j -0-^p^^(x',y';q)Fy,,ix',y';p-q), (6.13) 



*We note that this approximation can also be related to a nonperturbative 1/Nf 
expansion at next-to-leading order. 

^The relevant self-energies can be obtained with 

i5r2 ST2 SQab 

2« jT^ = 2z — — — = Sab ^ba = iV^ , 

0G4, oyAB 

and similarly for the fermionic self energies Ey = S 6ij : 

.sr2 



6D 
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Figure 1: Two- loop contribution to T2[G,D]. The solid and dashed lines 
represent the full fermion (D) and boson (Q) propagators, respectively. 



4' 



and for the fermion self-energies 



A^y{x^y';p) = -g'N, 



+p^y{x',y';q)F^ix',y';p-q) 
-0-^]^Fi^(x',y^;q)F^{x',y^;p-q) 



(6.14) 



(6.15) 



(6.16) 



Finally, we note that for the case of a non-vanishing scalar self-interaction, 
the only additional two-loop contribution gives rise to a local mass shift 
as described above. In particular, at this order there is no additional 
contribution to the direct scattering part, i.e. to the RHS of Eqs. fl6.10p 
and ()6.11|) relevant for thermalization. 



7 Initial conditions 

The time evolution for the fermions is described by first-order 
(integro-)differential equations for F and p: Eqs. ()5.5|) - ()5.8|) . As pointed 
out above, the initial fermion spectral function is completely determined 
by the equal-time anticommutation relation of fermionic field operators 
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(cf. Eq. (j3.13|) ). To uniquely specify the time evolution for F we have to 
set the initial conditions. The most general (Gaussian) initial conditions for 
F respecting spatial homogeneity, isotropy, parity, charge conjugation and 
chiral symmetry can be written as 



\t=t'=0 



1 
2 
0. 



- nl{p) 



(7.1) 
(7.2) 

Here nl{p) denotes the initial particle number distribution, whose values can 
range between and 1 (the definition of effective particle number distribution 
in terms of equal-time two-point functions is detailed in |Appendix A ) . At late 
times, when thermal equilibrium is approached, this will lead to a canonical 
description with zero chemical potential.^'' 

The evolution equations ()6.11|) and ()6.10|) for the scalar correlators are 
second-order in time and one needs to specify initial conditions for the 
propagators and their time derivatives. As for the fermions, the initial 
conditions for the scalar spectral function is completely specified by the field 
commutation relations (cf. Eq. ()6.9p ). For the scalar two-point function F^ 
we consider (cf. also Refs. 01221123) 



F^{t,t\p)\t=t'=o 

dtF^{t,t',p)\t=t>=o 
dtdt>F^{t,t',p)\t=tf=o 



eo(p) 
0, 



(7.3) 



eo(p) no{p) + 



with an initial particle number distribution no{p) and initial mode energy 
eo{p) (see Appendix A] ). 



It is instructive to consider for a moment the solution of the free field 
equations, which can be obtained from Eqs. ()5.5p - ()5.6|) and ()5.7|) - ()5.8|) by 
neglecting the memory integrals on their RHS. The solution of the fermionic 
free field equations with the above initial conditions reads 



Fv{t,t',p) 
F^it,t',p) 



nl{p)j cos[p{t - t')] , 
^ -^o(p)^ sm[p{t-t') 



^^In particular, for nonzero chemical potential or net charge density the BCS mechanism 
can lead to the condensation of Cooper pairs of fermions, which will be discussed elsewhere. 
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and for the spectral functions one obtains 



pv{t,t',p) = sm[p{t-t')] , p'y{t,t',p) = tcos[p{t-t')]. (7.4) 

One observes that each mode of the equal-time correlator Fy{t, t,p) is strictly 
conserved in the absence of the memory integrals. Since this correlator is 
directly related to particle number (see Appendix A), this means that the 



latter is conserved mode by mode in this approximation. Although this is 
expected in the free field limit, this is of course not the case in the fully 
interacting theory. We emphasize that such additional conservation laws do 
not only appear in the free field limit, but are a property of mean-field-type 
approximations, which include local corrections to the bare mass, but neglect 
the scattering contributions described by the memory integrals on the RHS 
of Eqs. ()5.5p - (j5.8|) . In these approximations, the existence of this infinite 
number of spurious conserved quantities prevents the system to approach 
the thermal equilibrium limit at late times. This aspect has been discussed 
in detail in the context of scalar field theories in Refs. It is therefore 

crucial to go beyond such "Gaussian" approximations in order to correctly 
describe in particular the late-time evolution of the system in the interacting 
theory. 



8 Numerical implementation 

We numerically solve the evolution equations (j5.5j) - (j5.8|) and (j6.1(J|) - (j6.1ip . 
together with the self-energies Eqs. ()6.13p - ()6.16|) . The structure of the 
fermionic equations is reminiscent of the form of classical canonical equations. 
In this analogy, Fv{t,t') plays the role of the canonical coordinate and 
Fy^t, t') is analogous to the canonical momentum. This suggests to discretize 
Fy{t,t') and pvit,t') at t - t' = 2nat (even) and F^{t,t') and p^(t,t') at 
t — t' = (2n + l)at (odd) time-like lattice sites with spacing at. This is a 
generalization of the "leap-frog" prescription for temporally inhomogeneous 
two-point functions. This implies in particular that the discretization in the 
time direction is coarser for the fermionic two-point functions than for the 
bosonic ones. This "leap-frog" prescription may be easily extended to the 
memory integrals on the RHS of Eqs. ()5.5|) - ()5.8p as well. 

We emphasize that the discretization does not suffer from the problem 
of so-called fermion doublers [25. The spatial doublers do not appear since 
fl5.5|) - ()5.8p are effectively second order in x-space. Writing the equations for 
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Fv(t,t',x) and pv{t,t',x) starting from (j5.5p - (|5.8|) one realizes that instead 
of first order spatial derivatives there is a Laplacian appearing. Hence we 
have the same Brillouin zone for the fermions and scalars. Moreover, time- 
like doublers are easily avoided by using a sufficiently small stepsize in time 
at/ as. 

The fact that Eqs. (j5.5j) - (j5.8|) and (j6.1(Jj) - (j6.11|) contain memory integrals 
makes numerical implementations expensive. Within a given numerical 
precision it is typically not necessary to keep all the past of the two-point 
functions in the memory. A single PHI desktop workstation with 2GB 
memory allows us to use a memory array with 470 timesteps (with 2 temporal 
dimensions: t and t'). We have checked for the presented runs that a 30% 
change in the memory interval length did not alter the results. For a typical 
run 1-2 CPU-days were necessary. 

The shown plots are calculated on a 470 x 470 x 32^ lattice. (The 
dimensions refer to the t and t' memory arrays and the momentum-space 
discretization, respectively.) By exploiting the spatial symmetries described 
in Section the memory need could be reduced by a factor of 30. We have 
checked that the infrared cutoff is well below any other mass scales and that 
the UV cutoff is greater than the mass scales at least by a factor of three. 

To extract physical quantities we follow the time evolution of the system 
for a given lattice cutoff up to late times and measure the renormalized 
scalar mass m which is then used to set the scale. In the evolution equations 
we analytically subtract only the respective quadratically divergent terms 
obtained from a standard perturbative analysis. We emphasize that for the 
results presented below we use the late-time thermal mass to set the scale, 
and not the vacuum mass for convenience. 

We made runs for a range of couplings g"^ = 0.49 - 1 which show very 
similar qualitative behavior. Below, we present plots corresponding to g = 1, 
for which the time needed to closely approach thermal equilibrium is the 
shortest. This allows us to obtain an accurate thermalization with the lowest 
numerical cost. 

9 Far-from-equilibrium dynamics 

In Fig. 121 we present the time evolution of the fermion equal-time two-point 
function Fv{t,t;p) for three momenta p. Results are given for two very 
different initial particle number distributions, which are displayed in the 
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5 10 15 20 25 30 50 100 500 

Figure 2: The time evolution of the fermion two-point function Fv{t, t; p) for 
three values of the momentum p, in units of the renormalized scalar thermal 
mass m. The evolution is shown for two very different initial conditions with 
the same initial energy density. One observes that the dynamics becomes 
rather quickly insensitive to the initial distributions displayed in the insets - 
much before the modes settle to their final values. The long-time behavior 
is shown on a logarithmic scale for t > 30m^^. 

insets (see also Eqs. (I7.1j) - (j7.3|) ). The (conserved) energy density is taken 
to be the same for both runs. In this case, since thermal equilibrium is 
uniquely specified by the value of the energy density, the correlator modes 
should approach universal values at late times if thermalization occurs. 

It is striking to observe from Fig. |21that after a comparably short time, 
much before the correlation modes reach their late-time values, the dynamics 
becomes rather insensitive to the details of the initial conditions: for a given 
momentum, the curves corresponding to the two different runs come close 
to each other rather quickly. During the slow subsequent evolution the 
system is still far away from equilibrium before the approach to the late- 
time values sets in. From both runs one observes that the characteristic time 
needed to effectively loose the information about the details of the initial 
conditions is much shorter than the time needed to approach the late-time 
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Figure 3: The same as in Fig. |21 but for the bosonic two-point function 
F^{t,t;p) for three different momenta. The initial particle number 
distributions for the two runs denote by "A" and "B" are displayed in the 
insets. For the employed parameters one observes that, in the scalar sector, 
the time needed to become effectively insensitive to the initial distributions 
is comparable to the time scale describing the approach to the universal 
late-time value. 



result. Moreover, the late-time values are found to be universal in the sense 
that the different runs agree with each other to very good precision. 

Fig. 121 shows the corresponding behavior of the scalar correlator modes 
F^{t,t;p). The respective initial particle number distributions in the scalar 
sector for the two runs are given in the inset. For the two different runs 
the modes having the same momenta approach each other rather slowly as 
compared to the fermionic sector. However, they reach their final values on 
a time scale which is comparable to that observed for the fermions in Fig |21 

To characterize these time scales in more detail, we consider in Fig. |3]the 
unequal-time two-point function Fv(t,t';p). As expected, if the system is to 
become insensitive to the details of the initial conditions, we observe that 
the correlation between some time t and another time t' is suppressed for 
sufficiently large t — t'. We note that the oscillation envelope of Fyit^t'^p) 
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Figure 4: The fermionic two-point function Fv{t,t';p) as a function oit — t' 
at late time t = 600 m^^. One observes a very good agreement with 
an exponential behavior. The rate is well described by the width of the 
corresponding spectral function in frequency space as shown in the inset (cf. 
the text for details). 



can be well described in terms of an exponential for sufficiently late times. 
The corresponding damping rate approaches a constant value. To estimate 
the asymptotic rate we show in Fig. |3]the unequal-time two-point function 
Fv{t,t';p) as a function of t — t' for the late time mt = 600. The fit to an 
exponential yields the damping rate 7j'^^™^^(p = 0.78m) = 0.03(1) m. We find 
a moderate momentum dependence of this rate with 'y^^'^^'^\o) = 0.067(l)m 
and 7f'"''^(0)>7r"''^(P>0). 

We emphasize that the rate 'y^^^^^^ (p) can be related to the width of the 
Fourier transform of the spectral function with respect to the time difference 
t — t'. In principle, the latter involves an integration over an infinite time 
interval. However, since we consider an initial-value problem for finite times 
we know pv{t,t',p) only on a finite interval. To overcome this problem, 
we fit the data for pv{t,t',p) by a 7-parameter formula that is capable to 
account for the observed oscillations and damping, but which is more general 
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than the usual 2-parameter Breit-Wigner formula. We perform the Fourier 
transformation on the extrapolated data. The resulting function py{uj,p) is 
displayed as a function of frequency uj in the inset of Fig. |3] One clearly 
observes a nonzero width of the spectral function, the value of which may 
be obtained from a fit to a Breit-Wigner formula. By doing so, we obtain 
a very good agreement of the damping rate inferred from the width of the 
spectral function on the one hand and from the linear fit on the log-plot for 
Fv(t,t';p) on the other hand. 

We can use the fermion damping rate to quantify the time scale 
characterizing the effective loss of the details of the initial conditions: 
Comparing with Fig. |21 we observe that the inverse fermion damping rate 
at p ~ ^1/71^^^™^^ = 15(l)m~^^ characterizes well the time for which the 

dynamics becomes rather insensitive to the initial distributions. In contrast, 
we find that this time scale does not characterize the late-time behavior. For 
the latter, we observe to very good approximation an exponential relaxation 
of each mode Fvit,t;p) to its universal late-time value. Carrying out the 
measurement for different modes, we observe that the corresponding rate is 
almost independent of momentum and is given by 1/7^'^™) = 95(5) 

The corresponding time scale 1/7!*^°'^™^ is therefore much larger than the 



A similar analysis can be performed for the scalar sector. Here we find 



at p ~ the corresponding values 1/7! ^ = 50(5) m ^ and I/7, 



90(5)m~^. One observes that although the damping rates for fermions and 
bosons are very different, the respective thermalization rates are rather 
similar. For the scalars the thermalization rate is larger than the damping 
rate, as is observed above for the fermions. However, the difference is 
much less pronounced for the scalars. Similar studies in 1 -|- 1 dimensional 
quantum jU] and classical |8J scalar theories typically find a substantial 
difference between damping and thermalization rates. However, this may be a 
consequence of the stringent phase-space restrictions and therefore particular 
to 1 + 1 dimensional systems. 

We finally note that approximate rates describing the early-time behavior 
may be determined by an exponential fit to the functions ^^(t, 0), F°(t, 0) 
and Fv{t,0) in a finite time interval. These rates depend on time and 
approach the late-time values given above. At early times we observe an 
approximate exponential damping with a rate about twice as big as the late- 
time value for the fermions, and about half the late-time rate for the scalars. 
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10 Approach to quantum thermal equiUbrium 

In the previous section, we have seen that the out-of-equihbrium evolution 
of the system leads to a universal late-time behavior, uniquely characterized 
by the initial energy density. We now analyze in detail if quantum thermal 
equilibrium, characterized by Bose-Einstein and Fermi-Dirac statistics, is 
approached^^ . In thermal equilibrium, the spectral function and the 
statistical two-point function are not independent of each other, but are 
related by the fluctuation-dissipation relation. The latter is an exact relation, 
which can be stated in 4-dimensional Fourier space as 122] 



for the scalar correlators, where n^^iu) = l/[exp(a;/T) — 1] denotes the Bose- 
Einstein distribution function. The frequency u is the Fourier conjugate of 
t — t' (in thermal equilibrium, time-translation invariance implies that two- 
point functions only depend on t — t'). The value of the temperature T is 
determined by the energy density of the system. The fluctuation-dissipation 
relation for the fermion correlators is given by the equivalent expression with 
the replacement {n^^{uj) + |) — ^ {n^^{uj) — i) in Eq. (jlO.lj) and the Fermi- 
Dirac distribution n^^ = l/[exp(u;/T) + 1]. The same relations hold as well 
for the statistical and spectral components of the self-energies in thermal 
equilibrium. They provide an unambiguous way to extract the distribution 
functions, which are specific for quantum thermal equilibrium. Out of 
equilibrium, the spectral and statistical two-point functions are completely 
independent in general. However, if thermal equilibrium is approached at 
late times, they become related by the fluctuation-dissipation relations. 

To extract the answer about the late-time distributions, without relying 
on any assumptions, we consider the statistical and spectral correlators in 
Wigner coordinates. For this we express F^{t,t']p) and p^{t,t';p), as well 
as the corresponding fermion two-point functions, in terms of the center 

"'^"'^We emphasize that thermal equiUbrium cannot be reached on a fundamental level from 
time-reversal invariant evolution equations at any finite time. The results demonstrate that 
thermal equilibrium can be approached very closely at sufficiently late time, without again 
deviating from it for practically accessible times. For a more detailed discussion of this 
aspect see Ref. 0. 
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Figure 5: The late-time ratio of the statistical two-point function and the 
spectral function in frequency space, both for fermions {Fy/py) and for 
scalars {Ftp/ptp). In thermal equilibrium the quotient corresponds to the 
Bose-Einstein (BE) distribution function for scalars and to the Fermi-Dirac 
(FD) distribution for fermions - independently of any assumption on a quasi- 
particle picture (see Eq. p0.1|) ). The BE/FD distributions are displayed by 
the continuous curves parametrized by the same temperature. The value 
of the latter, T = 0.94m, is actually not fitted but has been taken from the 
inverse slope of Fig. |H1 This shows the correspondence with the quasi-particle 
picture described in the text. 



coordinate X° = {t + t')/2 and the relative coordinate = t — t' and write 
p^{X\u,p) = / ds°e^-^°p^(X0 + s72,X0-s72;p), (10.2) 



-2X0 

and equivalently for the other correlators. Since we consider an initial-value 
problem, the time integral over s° is bounded by ±2X° (cf. also the detailed 
discussion in Ref . [22] ) ■ If thermal equilibrium is approached for sufficiently 
large , then the correlators do no longer depend on X° and a Fourier 
transform with respect iot — t' can be performed to very good approximation 
(cf. the discussion in Sect. Ej). The distribution functions may then be 
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Figure 6: The time-dependent fermion quasi-particle distribution nf(t,p) as 
a function of mode energy p for various times t. We have plotted the inverse 
slope function log(l/'n,/ — 1), which reduces to a straight line intersecting the 
origin when nf{t,p) approaches a Fermi-Dirac distribution. This plot shows 
the data for run "A" of Fig. El 

extracted from the quotient of the Wigner transformed two-point functions. 
In Fig. El we show the respective ratios at late times. Both functions are 
in good agreement with the equilibrium distributions n^^{uj) and n^^{uj), 
respectively (cf. Eq. ()10.H) ). We emphasize that the displayed continuous 
curves are no separate fits. They are the Bose-Einstein and Fermi-Dirac 
distribution functions parametrized by the same temperature. In particular, 
the value for the temperature is not fitted but has been extracted from the 
inverse slope parameter as explained below. 

We stress that the above procedure to extract the distribution functions is 
independent of any assumption about a quasi-particle picture. However, for 
many practical purposes it is very convenient to have an effective description 
of particle number and mode energy directly in real time - without the need 
of a Fourier transform. An efficient description is elaborated in jAppendix Aj 
The value for the temperature in the thermal equilibrium distributions of 
Fig. |31 has been actually measured based on this quasi-particle picture. In 
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Figure 7: The time-dependent boson quasi-particle distribution n(t,p) as a 
function of mode energy e(t,p) (see text) for various times. In this case the 
inverse slope function is log(l/ra + 1), which reduces to a straight hne in case 
of a Bose-Einstein distribution. This plot shows the data for run "A" of 
Fig. El 



[Appendix A| we define the effective particle number and energy to be used, 
both for the fermionic and bosonic fields. For scalar field theories, this quasi- 
particle picture has been successfully employed previously to investigate 
thermalization '^.'^l^. 

In Figs.|niand[Z|we show the effective quasi-particle number distributions 
defined as 

^-nj{t,p) = Fvit,t;p) (10.3) 

for fermions and 

l+n(t,p) = e{t,p)F^{t,t;p), (10.4) 

^'^^ V F^{t,t';p) )^^, ^ ' 

for scalars, where e{t,p) is the quasi-particle mode energy as discussed 
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Figure 8: Boson and fermion quasi-particle distributions at late times as 
a function of mode energy {n,p = n(t,p) and = e(t,p) for bosons and 
rip = nf{t,p) and €p = p for fermions). Both inverse slope functions employed 
in Figs, ini and [7| have been applied here in order to demonstrate that the 
fermions clearly do not follow a Bose-Einstein distribution, and the scalars 
cannot be characterized by a Fermi-Dirac distribution. In contrast, for the 
respective correct statistics the curves lie on top of each other, showing 
that both fermions and scalars are described by the same inverse slope 
parameter. The value of the latter is used in Fig. El as the temperature 
for the distributions n^^{u}) and n^^{uj). 



in [Appendix A The curves correspond to the initial conditions of run 
"A" shown in Figs. El and El One observes how the effective fermion 
and boson particle numbers change with time, the former approaching a 
Fermi-Dirac and the latter a Bose-Einstein distribution. To emphasize this 
point, we plot the corresponding "inverse slope functions" log(l/n/ — 1) and 
log(l/n+l), which reduce to straight lines for Fermi-Dirac and Bose-Einstein 
distributions respectively. The associated inverse slopes correspond to the 

^^Note that because of chiral symmetry there is no mass term present for the fermions 
and the corresponding quasi-particle mode energy is simply p. 
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temperature of the thermal equihbrium distributions. We see in Figs. E] and [7| 
that both inverse slope functions approach straight lines at late times. The 
associated temperatures for fermions and bosons are independent of time 
and agree very well with each other, as shown in Fig. |H1 For comparison, 
we display in Fig. |H1 the fermion inverse slope function evaluated with the 
bosonic effective particle number, and vice versa. This illustrates the degree 
of sensitivity of the inverse slope functions to the different statistics and in 
turn the degree of precision with which we are able to probe thermalization. 

11 Conclusions 

In this paper we have discussed the far-from-equilibrium dynamics and 
subsequent thermalization of a system of coupled fermionic and bosonic 
quantum fields. We solved the nonequilibrium dynamics beyond mean-field 
type approximations by calculating the complete lowest non-trivial order in 
a systematic coupling expansion of the 2PI effective action, which includes 
direct scattering as well as memory and off-shell effects. To our knowledge 
this is the first time that such a calculation is performed without further 
approximations. As a result, we show that, for various far-from-equilibrium 
initial conditions, the late-time behavior is universal and uniquely determined 
by the value of the initial energy density. Moreover, we are able to probe the 
approach to quantum thermal equilibrium, characterized by the emergence of 
Fermi-Dirac and Bose-Einstein distribution functions, with high accuracy. 
We emphasize that in the present calculation, besides the limitations of a 
coupling expansion, there is no other input than the dynamics dictated by the 
considered quantum field theory for given nonequilibrium initial conditions. 

This work can be extended in many directions. Most of the equations 
we derive are also valid in the phase with spontaneous symmetry breaking. 
Combined with earlier work on scalar theories jHl El I2n|; this provides a 
description of the nonequilibrium dynamics of the linear cx-model for QCD 
with two quark flavors. The model has served for many years as a valuable 
testing ground for ideas on the equilibrium phase structure of low energy 
QCD at nonzero temperature and density. With the present techniques a 
quantitative understanding of the out-of-equilibrium physics of this model is 
within reach. 
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Appendix A Effective particle number and energy 

A particle number can only be strictly defined in the presence of a conserved 
charge. However, for physical interpretation it is often convenient to define 
an effective particle number even when there is no conserved charge. In 
particular, besides a total particle number it is often useful to have a 
definition of an effective particle number per (momentum) mode. The latter 
is typically not conserved in an interacting theory, and in the context of 
thermalization one would like to find a time-dependent particle number 
distribution, which allows one to observe a Bose-Einstein or Fermi-Dirac 
distribution at sufficiently late times^^. Of course, the notion of an effective 
particle number is typically only meaningful if the relevant degrees of freedom 
or "quasi-particles" are weakly interacting. 

There are many different ways of how one can introduce the notion of 
an effective particle number (for a recent discussion see e.g. Ref. ,29^ and 
references therein). For example, one can define it as the average energy per 
mode divided by the energy of the corresponding mode. In an interacting 
theory, this procedure can be ambiguous because the expression of the total 
energy receives contributions from interactions and the average energy per 
mode is not uniquely defined. In this appendix, we discuss an elegant 
way to circumvent this difficulty. More importantly for our purposes, the 
procedure leads to a definition of an effective particle number density, which 
indeed allows one to directly observe the emergence of a Fermi-Dirac or 
Bose-Einstein distribution from the nonequilibrium dynamics for the theory 
considered in this paper. In particular, it can be explicitely shown that 
the effective particle numbers are always positive and, for the fermions, 
smaller than one (see [Appendix B] ). Applying the present construction to 

^■^We emphasize, however, that the definition of an effective particle number distribution 
is not necessary to analyse the approach to quantum thermal equilibrium, as is discussed 
in Sect.Uni 
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the case of neutral scalar fields, we recover the particle number definition 
used in previous studies to exhibit thermalization in purely scalar field 
theories 

Fermions 

We start by considering the case of charged fields and construct the 
effective particle number from the conserved current generated by the U{1) 
symmetry. The associated 4-current for each given flavor is ~ ip'^^^il). 
Fourier transformed with respect to spatial momenta, the expectation value 
of the latter can be written as J'^{t,p) = tr[7^D*^(t, t,p)], where the subscript 
stands for fermions (cf. Eq. ()3.1|) ). In terms of the equal-time statistical two 
point function, its temporal and spatial components read:^^ 

Jj(t,p) = 2[l-2F°(t,t;p)], 
Jf{t,p) = -4vFv{t,t-p). 

A nice property is that the above expressions in terms of the full correlators 
do not contain any explicit dependence on the interaction part. In order to 
obtain an effective particle number, we want to identify these expressions 
with the corresponding ones in a quasi-particle description with free-field 
expressions. These are given hy}^ 

j;(«^)(t,p) = 2[i +Qj{t,p)], 

Jr\t,p) = -2v[l~2Nf{t,p)], 

where Qf{t,p) = nj — ftf is the difference between particle and anti-particle 
effective number densities and Nf{t,p) = [uf + nf)/2 is their half-sum. The 
physical content of these expressions is simple: the temporal component J° 
directly represents the net-charge density per mode Qfit,p), whereas the 
spatial part J is the net current density per mode and is therefore sensitive 
to the sum of particle and anti-particle number densities. Identifying the 

^''As we shall see below for scalar fields, the expression for the effective particle number 
one obtains in this way can be directly applied to the case of neutral fields as well. 

^'''The constant factor in the temporal component conies from the fact that we define 
the current without the standard normal ordering. 

^^Here, we have explicitly used our assumptions of parity symmetry and rotational 
invariance (cf. Sect. |3J), which imply that the different spin states contribute the same, 
therefore the factor of two. 
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above expressions, we define 

^Qfit,p) = -F^{t,t;p), (A.l) 

^-Nf{t,p) = Fv{t,t;p). (A.2) 

Of course, tfiese expressions are only meaningful for the case that the 
interacting theory is well-described by a quasi-particle picture. The 
important point is that the above procedure allows one to construct an 
effective particle number density without knowing a priori which part of 
the interaction is to be considered as the "dressing" of the quasi-particles 
and which part describe their "residual" interactions. We note that the 
equal-time two-point functions on the RHS of these equations are real by 
definition (see Eq. ()5.3p ). Moreover, using the anticommutation relation for 
the fermion fields, it is shown in Appendix B| that these definitions always 



satisfy < Nf{t,p) < 1 and —1 < Qf{t,p) < 1. These properties are 
important for the above definitions to be physically meaningful. 

Introducing a non-vanishing net charge density per mode was useful for 
the above general construction. However, in the present paper, we consider 
only CP-invariant systems, which imply that the latter should vanish. 
Indeed, we see from Eqs. ()5.3p that the requirement of CP-invariance imply 
that our above definition of net charge density per mode vanishes identically 
for all times and for all modes. Therefore, the effective particle and anti- 
particle numbers are equal and we have 

Qf{t,p) = 0, 

--nf{t,p) = Fv{t,t;p). 

where nf{t,p) is the effective particle number. 
Bosons 

Following the same lines as above, let us consider for a moment the case of 
a single charged scalar field 0. The 4-current associated the corresponding 
f/(l) symmetry is ~ i[(f)^ {d^ (f)) — ((9^0^)0] and, as for the case of fermions, 
its expectation value has a simple expression in terms of the equal-time 
statistical two-point function of the charged field. In momentum space, it 
reads: 
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Jt,{t,p) = 2pF^{t,t;p), 

where, as before, the constant contribution comes from the fact that we define 
the current without the usual normal ordering. Here, the statistical two-point 
function for the charged scalar field is defined as the expectation value of the 
anticommutator of two fields operators: F^(t,t';p) = 0'''(t',p)]+). 
It has the symmetry property F^{t,t';p) = F^{t',t;p), so that the above 
expressions are real. The corresponding quasi-particle expressions read: 

j;^'^''\t,p) = Qb(t,p)-i, 
jr\t,p) = ^[i+2iv,(t,p)], 

where Qb{t,p) and Ni,{t,p) have the same meaning as before in terms of 
effective particle and anti-particle number densities and where €{t,p) is the 
quasi-particle energy. We therefore define: 

Qt,{t,p) = i{dt-dt')F^{t,t';p)\t'=t, (A.3) 
l + 2Nb{t,p) = 2e{t,p)F^{t,t;p). (A.4) 

Note that the RHS of both expressions arc real quantities, as they should if 
the LHS are to be interpreted as charge and quasi-particle number densities 
respectively. 

It remains to define the effective quasi-particle energy e{t,p)}'^ For this 
purpose, we use the free-field like expression for the average energy per mode, 
which in the case of a single charged scalar field can be written as 

dtdt'F4t,t';p)\t'=t + e''{t,p)F^{t,t;p) = e{t,p) [2n{t,p) + 1 

in terms of the statistical two-point function. Therefore, one obtains for the 
effective quasi-particle energy: 



"'^''Notc that this was not necessary for the fermionic case, because of our assumption 
of chiral symmetry, which prevents an effective mass term. In fact, one can repeat the 
following argument to obtain for the fermion effective quasi-particle energy ef{t,p) = p, 
as employed above. 
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It is straightforward to show that the combination of equal-time correlators 
appearing on the RHS of ()A.5|) is indeed positive. Moreover, using 
commutation relations for the scalar field, one can show that the effective 
particle number Nb{t,p), as given by Eqs. ()A.4|1 and ()A.5|) . is a positive 
quantity (see Appendix B] ). We note that for the case of CP-invariant 



systems, the effective charge density per mode ()A.3|) vanishes, as it should.^® 
When dealing with neutral scalar fields, as it is the case in the present 
paper, we can use the same formula derived above. Notice that this is 
consistent since in that case one has F^(t,t';p) = F^(t',t;p), which directly 
implies that Qb{t,p) = 0. Therefore, in the present paper, we use the 
definition 

^ + n{t,p) = e{t,p)F^{t,t;p) (A.6) 

together with ()A.5|) for the boson effective particle number n{t,p) and mode 
energy e{t,p) for each individual scalar species [221 E]- 



Appendix B 

Here we show that the combinations of equal-time two-point functions used 
in this paper to define effective particle number densities are always positive 
and, for the fermionic case, smaller than one. It is demonstrated to be a 
simple consequence of the (anti-)commutation relations of field operators. 

Fermions 

As a first exercise, we show that for any operators (f and if"^ satisfying the 
anticommutation relations 

=1 

and 

= \ = 0, 

one has 

< (^V) < 1, (B.l) 

where the brackets denote an average with respect to any density matrix. The 
left inequality above is trivially obtained for example by inserting a complete 
sum of states between the operators y?^ and if. One obtains a sum of positive 

^^This immediately follows from the behavior of the statistical propagator under CP- 
transformation: Fff,(t,t';p) > F^{t',t;p). 
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quantities which is of course positive. To show the second inequahty, we 
introduce the hermitian operator = (/'V- Using the anticommutation 
relations above, it is easy to see that A^^ = A^, from which it follows that 



An' 



( {n - {N))" ) = {N') - {N)' = {N) (l - (AT)) . 



It is clear that An^ > and we have just shown that (A^) > 0. One therefore 
conclude from the above equation that (A^) < 1, as announced. 

We now come to our effective particle number densities. Using Eqs. ()A.1|1 
and ()A.2|) and recalling that Qf = Uf — fif and Nf = {rif + nj)/2, we get for 
the effective particle and anti-particle number densities: 

l-n}{t,p) = Fy(t,t,p) + F^(t,t,p) = ^tr[(7o + t;-7)i^(t,t,p) 
i-n;(t,p) = Fy(t,t,p)-F°(t,t,p) = -^tr[(7o-t;-7)i^(t,t,p) 

In terms of the fermionic field operators ip{t,p) and which satisfy the 

anticommutation relations 



and 

one has 



Fit,t,p) = ^{mt,p),^{t,p)]_ 



Therefore, one can rewrite (from now on, we drop the explicit t and p 
dependence) : 



Now we introduce the operator 



70 + W-7 
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in term of which, 

a=l 

where we exphcitely wrote the sum over Dirac indices. For a given Dirac 
index, it is straightforward to check that 

= (l-7otT-7) =1 (B.3) 

and 

where we speciahzed to the Dirac basis to write the last equahty of 
Eq. ()B.3jl .^^ Using (jB.lj) for each individual a, we conclude from ()B.2|) that 

< nf{t,p) < 1 

for any time t and any momentum p. It is straightforward to repeat the 
above arguments to show that 

< < 1. 

Bosons 

For scalars, we first use Eqs. ()A.5|1 and ()A.5j) to rewrite 

l + n{t,p) = ^JF^{t,t■p) [dtdt'F4t,t';p)]t'=t. 

Recalling the definition of the statistical propagator in terms of field 
operators: 

F4t,t;p) = {4>\t,p)(P{t,p)), 

dtd,F^{t,t'-p)]t,=t = (nt(t,p)n(t,p)), 

where n(t,p) and n^(t,p) are the canonical momenta conjugate to (f)^{t,p) 
and (j){t,p) respectively (e.g. Il{t,p) = (t)\t,p)), we see that 

n(t,p)>o <=> (0t0)(ntn) > i, 

where we dropped in the notation the explicit time and momentum 
dependence of field operators. The second inequality is nothing but 
Heisenberg's uncertainty relation, a direct consequence of the equal-time 
canonical commutation relation s of field operators |;33j . 

""^^This is more convenient, but not necessary. It is simple to adapt the argument to any 
basis. 
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